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We present a cluster-based density-functional approach to model charge transport through molec- 
| ular and atomic contacts. The electronic structure of the contacts is determined in the framework 

£SJ ■ of density functional theory, and the parameters needed to describe transport are extracted from 

finite clusters. A similar procedure, restricted to nearest-neighbor interactions in the electrodes, has 
been presented by Damle et al. [Chem. Phys. 281, 171 (2002)]. Here, we show how to systemati- 
cally improve the description of the electrodes by extracting bulk parameters from sufficiently large 
metal clusters. In this way we avoid problems arising from the use of nonorthogonal basis functions. 
For demonstration we apply our method to electron transport through Au contacts with various 
CZ3 , atomic-chain configurations and to a single-atom contact of Al. 
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I. INTRODUCTION 
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Advances in the experimental techniques for manipulating and contacting atomic-sized objects have turned the 
, ^/ vision of molecular-scale electronic circuits into a realistic goali 2 -^^ This has intensified the interdisciplinary efforts 
to study charge transport in nanostructures. Ideally, the circuits would be constructed in a bottom- up approach 
with functional units and all the wiring on the molecular scale. To approach the goal, present-day experiments in 
the area of molecular electronics concentrate on measuring the current-voltage response of single molecules in contact 
to metallic electrodes. In these studies, also purely metallic atomic contacts and wires serve as important reference 
systems £ 

In order to support the experiments and to stimulate further technological advance, theoretical modeling of atomic- 
scale charge transport is needed. Here one faces the challenge to describe infinitely extended, low-symmetry quantum 
systems that may, in addition, be far from equilibrium and involve strong electronic correlations. While a complete 
theoretical understanding is still lacking, sophisticated ab-initio methods have been developed for approximate but 
parameter-free numerical simulations. In order to study the prototypical metal-molecule-metal systems or metallic 
• • . atomic contacts, many groups use density functional theory (DFT) combined with noncquilibrium Green's function 
(NEGF) techniques^^^^^^^^^^^L 2 ^ Some shortcomings related to the use of DFT in this context 
have been pointed out, and solutions are being sought i 23 ' 24 i 25 ' 26 On the other hand, from a practical point of view 
DFT presently appears to be one of the most useful ab-initio electronic structure methods, since studies of quantum 
transport require dealing with a large number of atoms. Furthermore the metal-molecule-metal contacts are hybrid 
systems, where the central regions fre quently behave rather insulator-like, while the electrodes are metallic. For more 
complete discussions we refer to Refs. 

The DFT approaches can mainly be divided into two types. In the first one, atomic-sized contacts are modelled by 
periodically repeated supcrcclls, and computer codes developed for solid-state calculations are employed ^ 15 ' 22 The 
use of periodic boundary conditions facilitates the electrode description. However, the conductance is determined for 
an array of parallel junctions and may be affected by artifical interactions between them. The second type is based on 
finite clusters and originates more from the chemistry communityi£ & 1 4 1 1 6 1 2 1 It has the advantage that genuinely single- 
atom or single-molecule contacts are described, and it makes possible investigations of molecules of large transverse 
extent. The drawback is typically the description of the electrode, since it is difficult to treat bulk properties based on 
finite clusters. Furthermore the coupling between the device region and the electrode can be complicated by finite-size 
and surface mismatch effects. 
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To arrive at an ab-initio DFT description it is necessary to treat the whole system consistently by using the same 
basis set and exchange-correlation functional everywhere. The problem of the cluster-based approaches regarding the 
electrodes is apparent, for example, from the work of Refs. H@(2l|, where the authors resort to a separate tight-binding 
parameterization obtained from the literature. 50 Damle et al. proposed to resolve this issue by extracting electrode 
parameters from finite clusters computed within DFT. 16-31 However, their treatment of the electrodes should be seen 
as a first approximation, since only couplings between nearest neighbor atoms were considered. Furthermore, they 
finally use energy-independent self-energies, which is well-justified only for electrode materials with a constant density 
of states (DOS) near the Fermi energy. 

In this work we present a cluster-based DFT approach for the atomistic description of quantum transport. We 
follow the ideas of Ref. l3l|, but place special emphasis on the treatment of the electrodes. In particular, we show 
that extracting electrode parameters from small metal clusters can lead to an unphysical behavior of the overlap 
of the nonorthogonal basis functions in fc-space. The description of the electrodes can be improved systematically 
by employing metal clusters of increasing size. Our implemenation is based on the quantum-chemistry package 
TURBOMOLE, which allows us to treat clusters of several hundred atoms. In this way we obtain an ab-initio 
formulation of quantum transport in atomic-sized contacts, where the whole system is treated on an equal footing. It 
has the advantage that we can employ high-quality quantum-chemical Gaussian basis sets, which are well-tested for 
isolated systems. 

The theoretical framework of our approach is presented in Sec. |TIJ Several technical details, related to the use of 
nonorthogonal basis functions and the electrode treatment can be found in Apps . 1X1 and [Bl To demonstrate the power 
of our methods we study in Sec. [TTT] the transport properties of atomic contacts of Au and Al. The choice of these 
materials is motivated by the fact that Au exhibits a rather energy-independent DOS near the Fermi energy, while 
Al does not. Furthermore, for these systems we can compare our results to the literature. We fi nd goo d agreement, 
and demonstrate the robustness of our results. Further applications have been presented in Refs. [32.33.34.35.3^. We 
summarize our results in Sec. IIVI 

II. THEORETICAL APPROACH 
A. Electronic structure 

Our ab-initio calculations are based on the implementation of DFT in TURBOMOLE 5.9.— By ab-initio we mean 
that the simulations require no system-specific parameters. Self-consistent DFT calculations of large systems are 
generally very time-consuming. TURBOMOLE, however, is specialized in handling such systems, and offers several 
possibilities to reduce the computational effort. Thus, one can exploit point group symmetries, including non-Abelian 
ones. In this way, the calculations speed up by a factor given by the order of the point group. Further options are 
the "resolution of the identity in J" (RI- J) 38 i 39 and the "multipole-accelerated resolution of the identity J" (MARI- 
J)^ which are both implemented in the ridft module of TURBOMOLE. The approximations help to reduce the 
effort to compute the Coulomb integrals J, which are particularly expensive to evaluate. With the help of the RI- J 
approximation, known also under the name "density fitting" , the four-center-two-electron integrals can be expressed 
as three-center-two-electron ones4^ Calculations are faster by a factor of 10-100 as compared to standard DFT, 
but equally accurate. The MARI-J technique concerns the Coulomb interactions between distant atoms. They are 
divided into a near-field and a far-field part, where the near field is treated with RI- J and the far field by a multipole 
approximation. Compared to RI-J, it can accelerate the calculations by another factor of 2-74^ DFT requires the 
choice of an exchange-correlation functional^ We select the generalized-gradient functional BP86, 43,44 which is known 
to yield good results for large metal clusters . 45 ' 46 ! 47 ' 48 To express the orbital wave functions, Gaussian basis sets of 
split-valence-plus-polarization (SVP) quality are used ) 38 ' 39 ' 49 which are the TURBOMOLE standard. Within the 
closed-shell formalism of DFT, total energies of all our clusters are converged to a precision better than 10 -6 a.u. In 
order to obtain ground-state structures, the total energy needs to be minimized with respect to the nuclear coordinates. 
We perform such geometry optimizations or "relaxations" until the maximum norm of the Cartesian gradient has 
fallen below 10~ 4 a.u. 



B. Transport formalism 

We compute transport properties of atomic-sized contacts using the Landauer-Biittiker theory and Green's functions 
expressed in a local, nonorthogonal basis i 50 ' 51 The local atomic basis allows us to partition the basis states \i, a) into 
left (L), central (C), and right (R) ones, according to a division of the contact geometry. In the basis states, a refers to 
the type of orbital at the position of atom i. For reasons of brevity we will frequently suppress the orbital index. The 
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Hamiltonian (or Kohn-Sham) matrix Hi a ^p = (i, a\ H \ j, f$), and analogously the overlap matrix = (i, a \ j, 

can thus be written in block form, 

( H LL H LC \ 
H=\ H CL H cc H CR . (1) 
\ H RC H RR J 

Both S and H are real-valued and are hence symmetric. In addition, we assume the C region to be large enough to 
have Slr = Hlb = 0. Within the Landauer-Biittiker theory^ the linear conductance can be expressed as 



G = j- I ,IE 



-(E), (2) 



where f(E,T) = {exp [{E — /x) / (k R T)] + 1} 1 is the Fermi function, and the chemical potential /i is approximately 
equal to the Fermi energy E Rl /i w Ep. Using the standard NEGF technique, the transmission function is given by 

t{E) = Tr [T L {E)G CC {E)T R (E)G CC {E)\ = Tr [S(E)t(E)] (3) 

with the transmission matrix 



t(E) = VW) G a cc (E) VTlJe). (4) 

Here we define the Green's functions 

G CC (E) = [ES CC - H cc - V L (E) - ^(E)}- 1 (5) 
and G cc = [G C c]\ tne self-energies 

J7 X (E) = (Hex ~ EScx) 9 r xx(E) (H XC - ES XC ) , (6) 

and the scattering-rate matrices 

r x (E)=-2hn[Z r x (E)}, (7) 

where g r X x 1S ^^e electrode Green's function for lead X = L,R. At low temperatures, the expression for the conduc- 
tance simplifies to 

G=^-r(E F ) = G J2 T n(E F ), (8) 

n 

with Go = 2e 2 /h the quantum of conductance and r„ the eigenvalues of tH. The latter are the transmission proba- 
bilities of the transmission channels nJ& Also other observables, such as the thermopower or the photoconductance, 
can be studied based on the knowledge of t(E)J£-£££2£1 

Information on the energetics of a system may help to identify conduction mechanisms. Such information can be 
extracted from the spectral density^ 

P (E) = ± [G r (E) - G a (E)] = -ilm [G r (E)} . (9) 

ITT 7T 

Using this, we define the local density of states (LDOS) at atom i and its decomposition into orbitals a via 

LDOS 4 (£) = ^LDOS M (£) (10) 

a 

LDOS ia (£) = (s^qccWS 1 ^) (11) 

In App.[X]we discuss the approximations involved in this definition. There we also consider further the issues related 
to the use of nonorthogonal basis sets to evaluate the single-particle Green's functions and the electric current. 
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Figure 1: Quantum transport scheme. The conduction properties of an atomic-sized contact (a) shall be studied. For this 
purpose it is divided into a C region and two semi-infinite L and R electrodes. Using a similar division as for the contact, 
information on the electronic structure of the C region {Sec, Hcc) as well as the CL and CR couplings (Scl, Hcl and Sen, 
Hcr) is extracted from the ECC (b). In order to obtain the self-energies E£ and EJj (a), the remaining task is to determine 
the electrode surface Green's functions g r LL and g^R- This procedure is described further below in the text. 



C. Implementation of the transport method 

1. Central system 

In order to determine the transmission function t(E) we need a practical scheme to obtain the necessary information 
on the electronic structure. In Fig. [1] we present our procedure. The goal is to describe the whole atomic-sized contact 
[Fig- HI a)] consistently, by treating the L, C, and R regions with the same basis set and exchange-correlation functional. 
We obtain the parameters Sec and Hcc as well as the couplings to the electrodes Sex and Hex with X = L,R from 
the extended central system (ECC) [Pig. Hfb)], in which we include large parts of the tips of the metallic electrodes. 
The division of the ECC into the L, C, and R regions is performed so that the C region is identical to that in Fig.[lja). 
The atoms in the L and R parts of the ECC [blue-shaded regions in Fig.QJa)] correspond to that part of the electrodes 
that is assumed to couple to C. The partitioning or division of the ECC is commonly made somewhere in the middle 
of the metal tips, and we will also refer to it as "cut". The electrodes [L and R regions in Fig. QJa)] are modeled as 
surfaces of semi- infinite crystals, described by the surface Green's functions g r X x- They are constructed from bulk 
parameters, extracted from large metal clusters. Further below we discuss in detail, how this is accomplished. 

In our approach, we assume the metal tips included in the ECC to be large enough to satisfy basically two criteria. 
First, all the charge transfer between the L and R electrodes and the C part of the contact should be accounted for. 
This ensures the proper alignment of the electronic levels in C with Ep. Second, most of the metal tips, especially 
the L and R parts of the ECC, should resemble bulk as closely as possible. In this way, we can evaluate the surface 
Green's functions by using bulk parameters of an infinitely extended crystal. Owing to surface effects caused by the 
finite size of the ECC, this can be satisfied only approximately. The mismatch between the parameters in the L and 
R regions of the contact and the ECC will thus lead to spurious scattering at the LC and CR interfaces. In principle 
this resistance can be eliminated systematically by including more atoms in the metal tips of the ECC. On the other 
hand, if the resistance in the C region is much larger than the spurious LC and CR interface resistances, they will 
have little influence on the results. 



2. Electrodes 



We extract bulk parameters describing perfect crystals from large metal clusters. The complete procedure, which 
aims at determining the surface Green's functions g r X x with X = L, R is summarized in Fig. [21 In this work we study 
exclusively electrode materials with an fee structure, of which Au and Al are examples. 

In a first step [Fig. Ufa)] we construct spherical metal clusters, henceforth called "spheres" . They are made up of 
atoms at positions {Rj\Rj = En=iin%i A \Rj \ < R s P here } with the standard fee primitive lattice vectors a n and the 
sphere radius R s P here _ Here, we will use the vector of integer indices j = {31,32, 3s) to characterize the atomic position 
Rj. We do not relax the spheres, but set the lattice constant ao to its experimental literature valued Increasing the 
radius R s P here should make the electronic structure in the center resemble that of a crystal. From the clusters we 
extract the overlap and Hamiltonian between the central atom at position and the neighboring ones at position j 
(including j = 0). This yields the matrix elements S^ h ^ e and Hj^ h ^ e , where a and f3 stand for the basis functions 
of the atoms at j and 0. For reasons of brevity, we will often suppress the orbital indices. Sjq ere and fj^P here are 
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sphere risphere - ^sphere rjsphere ri *.o(X) u(X) 

_'_ n d jO '"jO A jO' n jO y jO'jO 8xx 

symmetrize rotate construct 

(b) (c) (d) X=LR 

Figure 2: In order to obtain the electrode Green's functions g X x f° r l ea d X — L, R, we determine bulk parameters from large 
metal clusters. In a first step (a) we extract overlap and hopping elements, SJq ere , f{^P here j from the cluster's central atom to 
all its neighbors. They are (b) symmetrized by imposing the space-group of the electrode lattice. After (c) a rotation to adapt 
them to the orientation of the respective electrode, (d) gxx is constructed with the help of a decimation procedure. 

then matrices with appropriate dimensions. 

The bulk parameters Sjq, Hj Q need to satisfy the symmetries of the fee space-group [Fig. Hfb)]. While S^ here 
depends only on the relative position of atoms, surface effects due to the finite size of the fee clusters lead to deviations 
from the translational symmetry for H^ here . A rotation may still be necessary to arrive at parameters S^' , Hjq , 
which are adapted to the orientation of the electrode X = L,R [Fig. [DJc)]. Details on the symmetrization procedure 
and the transformation of crystal parameters under rotations are presented in App. |5J The parameters S^q \ 
are finally employed to construct the semi-infinite crystals and to obtain the surface Green's function g xx [Fig. Hid)] . 
Due to the finite range of the couplings Sex, Hex, we need to determine g r X x f° r the first few surface layers only 
[blue-shaded regions in Fig. HJa)]. We compute these with the help of the decimation technique of Ref. [5?], which we 
have generalized to deal with the nonorthogonal basis sets^ The parameters Sjo, Hj can be computed once for a 
given metal and can then be used in transport calculations with electrodes of various spatial orientations. 

For Au (a = 4.08 A) we have analyzed spheres ranging between 13 and 429 atoms, while for Al (cto = 4.05 A) 
they vary between 13 and 555 atoms. Since we want to describe bulk, the parameters extracted from the largest 
clusters will obviously provide the best description. There is, however, an additional criterion, which necessitates the 
use of large metal clusters for a reliable description of the electrodes. As discussed in App. IB 11 it is based on the 
positive-semidefiniteness of the bulk overlap matrix. We find a strong violation of this criterion, if the extraction of 
parameters is performed such that only the couplings of the cental atom to its nearest neighbors are considered. As a 
further demonstration of the quality of our description, we show in App. lB"2l the convergence of the DOS with respect 
to R s P here . 

For the transport calculations we need a value for the Fermi energy. The biggest Au and Al spheres computed, 
AU429 and AI555, respectively, are very metallic. They exhibit differences between the highest occupied molecular 
orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of less than 0.06 eV. Therefore we set Ep 
midway between these energies. In this way we obtain Ep = —5.0 eV for Au and Ep = —4.3 eV for Al. The values 
will be used in all the results below. Notice that the negative values of Ep agree well with experimental work functions 
of 5.31 to 5.47 eV for Au and 4.06 to 4.26 eV for Al. 59 

III. METALLIC ATOMIC CONTACTS 

In this chapter we explore the conduction properties of metallic atomic contacts of Au and Al. These systems, 
in particular atomic-sized Au contacts, have been studied in detail both experimentally and theoretically, and can 
therefore be used to test our method. We start by discussing the transport properties of the Au contacts, consisting 
of a four-atom chain, a three-atom chain, and a two-atom chain or "dimer". Since Al does not form such chains, we 
consider only a single- atom contact. For all systems we analyze the transmission, its channel decomposition, and, 
in order to obtain knowledge about the conduction mechanism, the LDOS for atoms in the narrowest part of the 
contact. Moreover, we demonstrate the robustness of our transmission curves with respect to different partitionings 
of the large ECCs. 

A. Gold contacts 

Let us first discuss the electronic structure of Au, where we display the DOS in Fig. [3] The Fermi energy at 
Ep = —5.0 eV is located in a fairly structureless, flat region somewhat above the d band. Based on the electronic 
configuration [Xe].4/ 14 .5d 10 .6s 1 of the atom, one might have expected a strong contribution only from the s orbitals 
at Ep. But, as is visible from Fig. \Mb), s, p, and d states all yield comparable contributions^ This signifies that 
valence orbitals hybridize strongly in the metal. 
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Figure 3: DOS for Au. (a) DOS resolved into s, p, and d contributions and (b) DOS resolved into all individual orbital 
components. The dashed vertical line indicates Ef = —5.0 eV. 



When an atomic contact of Au is in a dimer or atomic chain configuration, a conductance of around IGq is expected 
from experimental measurement o 60 ' 61 ! 62 ' 63 ^ 64 as well as from theoretical studies.— i 65 i 66 i 67 i 68 The analysis shows that 
this value of the conductance is due to a single, almost fully transparent transmission channel. It arises dominantly 
from the s orbitals of the noble metal Au, since the electronic structure in the narrowest part of the contact resembles 
more the electronic configuration of the atom i 64 i 65 



1. Determination of contact geometries 

Despite the consensus that the conductance of atomic chains of Au is around IGo, the precise atomic positions 
play an important rolej 68 ' 69 Therefore it is necessary to construct reference geometries that have been studied with 
a well-established transport method. We choose to compare to results obtained with TRANSIESTAfI The ECCs 
investigated are shown in Fig. [4j The four-atom Au chain with electrodes oriented in the [100] direction, called 
Aul00c4, corresponds to a contact geometry examined in Ref. l70l [see Fig. 1(b) therein]. The three-atom Au chain, 
Aulllc3, is similar to a configuration in Fig. 9(d) of Ref. .7:. In addition, we study a Au dimer contact, Aulllc2, where 
a two-atom chain is forming the narrowest part. In contrast to Aul00c4, for the latter two contacts the electrodes are 
along the [111] direction. 

Let us briefly explain, how we determine these geometries (Fig. 2]). For Aul00c4 we construct two ideal, atomically 
sharp Au [100] pyramids, with two atoms in between. The pyramids end with the layer consisting of 25 atoms. The 
distance between the layers containing four atoms is set to 12.68 A [Fig. |Ua)], as in Ref. |7(J. Next we relax the four 
chain atoms without imposing symmetries, keeping all other atoms fixed. After geometry optimization, we find that 




(a)Aul00c4 (b)Aulllc3 (c)Aulllc2 



Figure 4: ECCs for Au. (a) Aul00c4 is a four-atom Au chain, (b) Aulllc3 is a three-atom Au chain, and (c) Aulllc2 is a 
two-atom Au chain or dimer contact. For (a) electrodes are oriented in the [100] direction, while for (b) and (c) this is the [111] 
direction. The main crystallographic direction is assumed to be parallel to the z axis. Indicated are also the most important 
bond distances together with information on the construction of the ECCs. 
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Figure 5: Aul00c4. (a) ECC with two different partitionings into L, C, and R regions, cuts 1 and 2, and (b) the transmission 
as a function of the energy for these cuts. For cut 2 (c) transmission resolved into its transmission channels and (d) LDOS of 
the chain atom indicated in (a). 



the configuration agrees well with symmetry D^, We add two more Au layers with 16 and 9 atoms on each side, where 
the ECC now consists of 162 atoms, and perform a final DFT calculation, exploiting the symmetry D^h- As compared 
to Ref . l70l. all bond distances indicated in Fig.[4ja) agree to within 0.01 A, except for the distance between the central 
chain atoms, where our distance is shorter by 0.07 A. For Aulllc3 we proceed similarly to Aul00c4 [Fig. HJb)]. We 
start with two perfect Au [111] pyramids, set the distance between the Au layers with 3 atoms to 9.91 Aj£ and cut 
the pyramid off at the layers containing 10 atoms. Then we add one atom in the middle, relax the three chain atoms, 
add two layers on each side with 12 and 6 atoms, and perform a calculation in symmetry D^d- Aulllc3 consists of 77 
atoms in total. Our bond distances agree with those in Fig. 9(d) of Ref. to within 0.02 A. For Aulllc2 we include 
also the first Au layer in the geometry optimization process. The distance between the fixed layers with 6 atoms is 
12.12 A. Otherwise the steps are the same as for Aulllc3. The ECC is computed in symmetry D$d and consists of 76 
atoms. In the parts excluded from the geometry optimization, atoms are all positioned on the bulk fee lattice, where 
we set the lattice constant to the experimental value of 4.08 A, which corresponds to a nearest-neighbor distance of 
2.88 A. In each ECC the main crystallographic direction is aligned with the z axis, which is the transport direction. 



2. Four-atom gold chain 

Let us now study the conduction properties for the contact Aul00c4 [Fig. OJa)]. There are different possibilities to 
partition the ECC into the i, C, and R regions. The cuts should be done so that L and R are unconnected (Slr = 
and Hlr = 0, see Sec. Ill Bp . Hence the C region must be long enough. In order to describe well the coupling to the 
electrode surface (Fig. [1}, it is furthermore necessary to have sufficiently many layers in the L and R regions. We 
observe that at least two layers are needed to obtain reasonable transmission curves. For the two different cuts of 
Fig-EJa) t(E) is plotted in Fig.[5jb). In both cases it is found to be almost identical, indicating a sufficient robustness 
of our method. The transmissions at the Fermi energy are r(Ep) — 0.93 and 0.98 for cuts 1 and 2, respectively. These 
values correspond well to the result r(Ep) — 0.99 of Ref. [70. 

For cut 2 it is visible in Fig. [5fc) that the transmission at Ep is dominated by a single channel, in good agreement 
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Figure 6: Aulllc3. (a) ECC with two different partitionings into L, C, and R regions, cuts 1 and 2, and (b) the transmission 
as a function of the energy for these cuts. For cut 1 (c) transmission resolved into its transmission channels and (d) LDOS of 
the central chain atom indicated in the ECC. 

with experimental observations^ and previous theoretical studieSfi^ In general, the electronic structure at the 
narrowest part should have the most decisive influence on the conductance of an atomic contact. Therefore we plot in 
Fig. G3^d) the LDOS of the atom indicated by the arrow in Fig. 03a), resolved in its individual orbital contributions. 
Compared to the bulk DOS of Fig. 03 it is dominated by s at Ep, where the contributions of all other orbitals than 
s and p z are suppressed. These two orbitals will form the almost fully transparent transmission channel, which is 
radially symmetric with respect to the z axis. 

3. Three-atom gold chain 

Exactly the same analysis will now be carried out for the contact Aulllc3. In Fig. [B]the geometry of the ECC, the 
transmission for different partitionings, the transmission channels, and the LDOS of the central chain atom are shown. 
As for Aul00c4, we observe that the different cuts yield very similar transmission curves [Fig. OUb)]. Furthermore 
all the basic features in t(E) are the same as in the TRANSIESTA calculation [see Fig. 11(d) of Ref.0]. Above the 
d band, which exhibits a very narrow and high final peak, there is a dip in t{E) in both cases. The transmission 
recovers, however, and a flat region with a value of around one is visible. At the Fermi energy cuts 1 and 2 yield 
r(Ep) = 0.96 and 0.99, respectively. This is in reasonable agreement with r(Ep) — 0.94 in Ref. 0, considering the 
differences in the electrode geometry, basis set, and exchange-correlation functional. We observe from Fig. [6jc) for 
cut 1 that the transmission at Ep is dominated by a single transmission channel, and the LDOS indicates a dominant 
contribution of s orbitals [Fig.^d)]. In addition, the peak structures in t(E) for the d states correspond well to peaks 
in the LDOS. This observation can also be made in Figs. 03c) and03!d) for Aul00c4. 
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Figure 7: Aulllc2. (a) ECC, and for the cut and atom indicated (b) the transmission resolved into its transmission channels 
and (c) the LDOS with its individual orbital contributions, respectively. 



4- Two-atom gold chain 



The transmission and LDOS resolved into transmission channels and orbital components, respectively, are shown 
in Fig. [7] for the dimer contact Aulllc2. As for Aul00c4 and Aulllc3 we observe a single dominant channel at Ep, 
and r(Ep) = 0.96. The finding of such a dominant channel for chains of two or more atoms is in good agreement with 
our analysis of less symmetric contacts, which were based on a combination of a tight-binding model and classical 
molecular dynamics simulations^ However, that t(E) increases partly even above one in the vicinity of Ep , signals 
that the influence of other channels is increased as compared to Aul00c4 and Aulllc3. Indeed, the LDOS of the atom 
in the narrowest part of the constriction [Figs. [7(a) and[7lc)] shows in particular increased p x and p v contributions. 
Also, the d states exhibit a less pronounced peak structure than was visible in Figs. Od) and [Hid). This is due to the 
higher coordination number of the atom and the enhanced coupling to the electrodes. 



B. Aluminum contacts 



As is visible from the bulk DOS in Fig. [51 the electronic structure of Al differs substantially from that for Au. While 
the latter is a noble metal with an s valence, the Al atom has the electronic configuration [NeJ.Ss 2 ^^ 1 with an open 
p shell, and the metal is hence considered sp-valent. The strong contribution of s and p states is also observed in the 
DOS, where d states play only a minor role. As compared to Au, the DOS exhibits a noticeable energy dependence 
around Ep. 

For Al we study an ideal fee [111] pyramid, consisting of 251 atoms [Fig. [H](a)] , henceforth referred to as Allllcl. 
Ideal means that the atoms are positioned on an fee lattice with the experimental lattice constant ag = 4.05 A. We 
have already reported results for Al dimer contacts in Ref. 1341 . 




Figure 8: DOS for Al. (a) DOS resolved into s, p, and d contributions and (b) DOS resolved into all individual orbital 
components. The dashed vertical line indicates Ef = —4.3 eV. 
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Figure 9: Allllcl. (a) ECC with three different partitionings into L, C, and R regions, cuts 1, 2, and 3, and (b) the transmission 
as a function of the energy for these cuts. For cut 2 (c) the transmission resolved into its transmission channels and (d) the 
LDOS of the central atom, indicated by the arrow in (a). The nearest-neighbor distance shown in the ECC is identical for all 
atoms. 



1. Aluminum single-atom contact 



In Fig. [9] the transmission is displayed for three different partitionings of the ECC Allllcl. Also shown are the 
transmission channels and the LDOS of the atom in the narrowest part of the contact for a selected cut. For energies 
below —6 eV, there are practically no differences visible between the curves for the three different partitionings. 
Nevertheless, some deviations arise at Ep, and we obtain r(Ep) — 2.36 (cut 1), 1.88 (cut 2), and 2.23 (cut 3). Similar 
to Au, we attribute these 20% relative variations to spurious scattering at the LC and CR interfaces. Our values 
for r[Ep) of around two agrees nicely with those reported for single-atom contacts in Ref. [HI- Compared to Au, the 
transmission-channel structure has ch ang ed in an obvious way. There are three channels at Ep, which is in line with 
experimental observations of Refs. [r34ll7ll . Due to the D^d symmetry of the ECC, transmission-channel degeneracies 
arise, where in particular r 2 — r 3 . As is visible from the LDOS, these additional channel contributions mainly stem 
from the p x and p y orbitals, while s and p z are forming the nondegenerate Tli 64 i 65 i 72 



IV. CONCLUSIONS 



We have developed a cluster-based method to study the charge transport properties of molecular and atomic 
contacts. We treat the electronic structure at the level of DFT, and describe transport in terms of the Landauer 
formalism expressed with standard Green's function techniques. Special emphasis is placed on the modeling of the 
electrodes and the construction of the associated bulk parameters from spherical metal clusters. We showed that 
these clusters need to be sufficiently large to produce reliable bulk parameters, where a criterion for the extent of 
the spherical clusters is set by the overlap of the nonorthogonal basis functions. In our studies we crucially rely on 
the accurate and efficient quantum-chemical treatment of systems consisting of several hundred atoms, made possible 
by use of the quantum chemistry package TURBOMOLE. Compared to supercell approached, our method has the 
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advantage that we genuinely describe single-atom or single-molecule contacts. 

As an application of our method we analyzed Au and Al atomic contacts. Studying a four-, a three-, and a two- 
atom chain with varied electrode lattice orientations for Au, we found a conductance close to IGo, carried by a single 
transmission channel. Next we investigated an ideal Al single-atom contact, and found three transmission channels 
to contribute significantly to the conductance of around 2Gq. These results are in good agreement with previous 
experimental and theoretical investigations. Both for Au and Al we demonstrated the robustness of our transmission 
curves with respect to partitionings of the contact systems. The results illustrate the applicability of our method to 
various electrode materials. 

Beside the metallic atomic contacts examined here, the presented method has been applied in the field of 
molecular electronics. Studies include the dc conduction properties of dithiolated-oligophenylene and diamino- 
alkanc junction a 32 i 33 ' 35 i 36 as well as oxygen adsorbates in Al contacts^ In addition, the thermopower— and 
photoconductance 3 - 3 . of molecular junctions has been investigated in this way. Our studies demonstrate the value 
of parameter-free modeling for understanding transport at the molecular and atomic scale. 
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Appendix A: NONORTHOGONAL BASIS SETS 

For practical reasons one often employs nonorthogonal basis sets in quantum-chemical calculations, consisting 
for example of a finite set of Gaussian functions. The electronic structure is described in the spirit of the linear 
combination of atomic orbitals (LCAO) ) 41 ' 73 ? 74 and this is also how TURBOMOLE is implemented. While it is in 
principle always possible to transform to an orthogonal basis, it may be more convenient to work directly with the 
nonorthogonal states. 

A concise mathematical description using nonorthogonal basis states can be formulated in terms of tensors. The 
formalism is presented in a fairly general form in Ref. l75l . where also the modifications of second quantization are 
addressed. Below we discuss some of the subtleties related to the use of nonorthogonal basis functions that are 
important for our method^ Since the basis functions are real-valued in our case, the full complexity of the tensor 
formalism is not neededi 76 i 77 Furthermore we use a simplified notation, where all tensor indices appear as subscripts 
of matrices. 



1. Current formula for nonorthogonal, local basis sets 

The most important quantity for transport calculations is the electric current. In the NEGF formalism, its deter- 
mination requires a separation of the contact into subsystems similar to Fig. [T](a) i 52 i 78 i 79 However, due to the overlap 
of the basis functions in a nonorthogonal basis, the charges of the subsystems are not well defined. Different ways of 
determining them exist, e.g. the Mulliken or Lowdin population analysis^ 3 - Despite these additional complications, 
the Landauer formula [Eq. ^] can be derived in a similar fashion as for an orthogonal basis. Recent discussions of 
the derivation can be found in Refs. [5lll80l 

2. Single-particle Green's functions 

Consider the single-particle Hamiltonian H describing the entire system. The retarded Green's operator is defined 
as G r {E) = [(E + i0 + ) 1 — H] 1 . Now consider the local, nonorthogonal basis \i) with the (covariant) matrix elements 
of the overlap Sij — (i \j) and the Hamiltonian Hij = (i\ H |j) . 76 i 77 Compared to Sec. Ill Bl the index i, used throughout 
this appendix, is a collective index, denoting both the position at which the basis state is centered and its type. The 
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components of the retarded Green's function, defined by G r = J2i i \i)GL(j\r^ satisfy the equation^ 



^2[(E + iO+)S jk -H jk ] G r kl (E) 



Ojl. 



(Al) 



The Green's function G cc is defined as G\j restricted to the central region C. It can be calculated according 
to Eq. ([5|). Due to the nonorthogonal basis, the perturbation that couples C to the lead X = L,R and enters the 
self-energy [Eq. ([6|], is given by Hex — EScx and thus includes also an overlap contribution. It is interesting to 
observe that as E — > oo the self-energies and the Green's function behave as 



X r x (E) E_ 



oo EScx (Sxx) 1 Sxc 



G 



CO 



E -> oo E- 1 (S~ 



'cc 



(A2) 
(A3) 



with 



cc 



Sec + ^ Sex {Sxx) 1 Sxc 



(A4) 



X=L,R 

Thus also the inverse overlap matrix of C is "renormalized" due to the coupling to the leads. 



3. Local density of states 



Using a set of orthonormal energy eigenstates \p) that satisfy H \p) = e M |/x), we obtain the decomposition p^ u (E) = 
(p\p(E)\v) = 6(E — e^S^v of the spectral density defined by Eq. ©. Clearly p fiv (E) fulfills the normalization 

/oo 
dEp llv {E) = 8 llv . (A5) 
- OO 

If, instead, we consider the components denned by pij(E) — — Im /V, where G\AE) is given by Eq. (|A1[) . we 

find 

/oo 
dE Pij (E)= (A6) 
-oo 

The normalization of Eq. (|A5|) can be recovered by performing a Lowdin orthogonalization of the basis 

/oo 
dE{S l / 2 ) %k p kl {E){S l ' 2 ) l0 = Sa (A7) 
-oo 

Let us analyze the LDOS of the central region pcc(E) = — Im \G r cc {E)] /it, which is Pij{E) restricted to C. 
Analogously to Eq. (|A7[) . we have defined the LDOS at atom i and its decomposition into orbitals a in Eqs. (fT0|) 
and (TTT|) . Since pcc(E) is a positive-semidefinite matrix, it is easy to show that LDOSi a (E) is positive for all 
E. However, the normalization dELDOSi a (E) = 1 is only approximately fulfilled. This could be corrected 

by multiplying in Eq. (fTT|) with (S^ 1 )^ 2 [Eq. (|A4I) ] instead of S c ^. But since the self-energy contributions 

J2x=l R^cx (Sxx) 1 Sxc constitute only a surface correction, their neglect may be justified for atoms in the 
middle of C. 



Appendix B: DESCRIPTION OF ELECTRODES 
1. Size requirement for the cluster construction 

How large do the spherical metal clusters, involved in the construction in Fig. [2 need to be for a convergence of 
the bulk parameters? Since the matrix elements of the Hamiltonian and the overlap decay similarly with increasing 
interatomic distance, we can concentrate on the overlap. For it a rather well-defined criterion can be found: The 
clusters should so large that the extracted bulk overlap matrix is positive-semidefinite. 
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Figure 10: s-orbital-chain model, (a) The overlap Sj Sl a s of an atom with its neighbors at positions Rj jX = jiao with ji = 
0, ±1, . . .. (b) Overlap S ss (k x ) after a discrete Fourier transformation. In both cases the solid line is for a large and the dashed 
line for a small cluster. 



We define states fc, aj — J2j e lk ' Rj \ j,a) in fc-space. Since Siajp = (i,a \ j,/3) is a positive-semidefinite matrix^ 
the same is true for the overlap in fc-space 

S a p{k, fc') = (fc, a |fc', P) =J2 e'^'Slcmpe^ (Bl) 

l,rn 

= NS^pS a/ 3(k), 

where we used that Si atm f3 — 5 , (;- m )Q,0/3- I n the expression, N is the number of atoms in the crystal and 

s a p(k)=j2 e ~ in '* iS J<*,w- ( B2 ) 

In order to study the positive-semidefiniteness of S a p(k,k') it is hence sufficient to investigate the behavior of 
Sapik). To do so for a complex quantum-chemistry basis set, we define the positivc-dcfinitcncss measure 

^ R sphere^ = min ( 5 (fc)). ( B3 ) 
k 

In this expression S(k) is the smallest eigenvalue of the matrix S a /s(k), where iS^^fc) is constructed from the crystal 
parameters extracted from a cluster with radius R s P here [S a fi(k) — J2r ••In -\<R>phere e ~ lk ' Ri Sja,op]- In the discrete 
Fourier transformations we assume periodic boundary conditions with a finite periodicity length along the standard 
primitive lattice vectors . 56 ' 58 R s P here must be chosen large enough for £ to be positive or, if £ remains negative, it 
must at least be sufficiently small in absolute value. 

Let us first illustrate the behavior of £ at the example of an s-orbital model. Gaussian s functions are described by 



Ur f )=(- S \ e-°W (B4) 



3/4 

IT J 

with an exponent a, characterizing the radial decay. Hence, the overlap between two atoms 

S js , 0s = [d 3 rMr~ R 3 )Mr) = e~"^ /2 (B5) 



decays with their distance Rj = \Rj\ like a Gaussian function. We consider an infinitely extended chain with atoms 

at equally spaced positions along the x-axis (Rj = jiaoe x ). The overlap from a selected atom to its neighbors drops 
off exponentially as shown in Fig. I10( a). The Fourier transformation will again result in a Gaussian with purely 
positive values S ss (k x ) [Fig. HOf b)]. If, however, overlap matrix elements are taken into account only up to a certain 
maximum value \Rj\ < R s P here j as in a finite cluster, a rough sin(k x )/fca;-behavior results, where S ss (k x ) becomes 
negative at certain fc- values. Upon an increase of R s P here i S ss (k x ) will evolve into a Gaussian function and £ will thus 
approach zero from below. The negative tails of S ss (k x ) are unphysical, and our observation implies that the clusters 
used to extract bulk parameters (Fig. ^ need to be of a sufficiently large radius R s P here i i n order to obtain a reliable 
description of a crystal. Obviously, the magnitude of R s P here depends on the basis set chosen. 



14 




12345 012345 



R spher 7a R sphere /a 

Figure 11: The positive-definiteness measure £ for Au and Al as a function of R s P here . Beside the SVP basis set, the behavior 
of £ is also shown for LANL2DZ, used in Refs.[]Jl3l S(k) [Eq. ((B2j] is evaluated at 32 3 fc-points. The radius R"P hETE has been 
scaled with the respective lattice constants (ao = 4.08 A for Au and ao = 4.05 A for Al). 



In Fig. [TT] we plot the behavior of £ as a function of R s P here for Au and Al. Beside the results for the SVP basis 
set, we display £ for Au also for the basis set LANL2DZ, used in Refs. H3l3ll It is visible that £ is positive for a single 
atom (R s P here = o), but negative for small spheres. With increasing R s P here j £ approaches from below similar to the 
s-orbital model. We find that the elimination of diffuse functions reduces the radius R s P here for £ to become positive 
or negligibly small. For practical reasons it may happen that R s P here cannot be chosen large enough to fulfill the 
positive-semidcfiniteness criterion £ > 0. In such a case negative eigenvalues of S a p(k) can lead to negative eigenvalues 
of the scattering-rate matrices Tx [Eq. ©], since pxx = —lm [g r X x\ / n ma Y 110 longer be positive-semidefinite (see 
also the discussion in Sec. I A 3[) . 



2. Bulk densities of states 

The DOS can be used as another measure for the convergence to a solid-state description. With a fc-space Hamil- 
tonian in an orthogonal basis set H orth (k), it is given as 

DOS(£) - VDOS Q (£) = --Tr a flm [g£$£(£)11 , (B6) 

1 — ' 7T 



where a runs over all basis functions on a bulk atom and G°^'\E) = J BZ d 3 kG orth ' r (k, E)/V B z with G orth ' r (k, E) = 
_ -I -l 

El - H orth (k) and with the volume V BZ of the first Brillouin zone BZ. The orthogonal Hamiltonian can be 
obtained in several ways. Two possible choices are (i) to Fourier transform Sjo and Hjo and perform a Lowdin 
orthogonalization in /c-space H orth {k) = S-y 2 (k)H(k)S-^ 2 (k) or (ii) to construct H$ th , which involves the Lowdin 
transformation H s P here - orth = ^ S ^/ 2 ) sphere H s P here (5-1/2) "* epe in real spac6j extraction of H sphere,orth &nd thc 
imposing of the fee space group, and to carry out the Fourier transformation only thereafter. For parameters extracted 
from large enough clusters we observe the equivalence of the DOS construction with respect to the two different 
orthogonal Hamiltonians H orth (k). If £ remains (slightly) negative due to a too small R s P here ; then the construction 
of the DOS from H°Q th [procedure (ii)] is of a higher quality than that resulting from the Lowdin orthogonalization 
in fc-space [procedure (i)]. 

In Fig. [12] we show the DOS as constructed via procedure (ii) with parameters extracted from different Au and Al 
spheres with 141 to 555 atoms. We observe that the DOS seems well converged with respect to R s P here both for Au 
and Al for the largest spherical clusters AU429 and AI555. 



3. Transformation of electrode parameters under rotations 



We assume that two coordinate systems are connected by the rotation g, where f* = gf. The transformation 
properties of the electrode parameters 



Y ja>ofi = (j,a|r|0,/3) = / dV«(i*--R*)iW/j(»D 
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with Y = S, H are determined by those of the basis functions (r | j, a) — 4> a {r — Rj)- The Gaussian basis functions 
used by TURBOMOLE are characterized by the angular momentum I and the multiplicity v — 1, . . . , 21 + 1, and a 
is a collective index for both. The rotated basis functions of angular momentum / can be expressed as 81 

22 + 1 



with the representation D l ^(g) of the rotation g. Using Y'(r) 
parameters of the two coordinate systems are related by 



Y(g 1 r), it can be shown that the electrode 



Y- - 



(B7) 



where D(g) is the representation of g in the employed basis set. By knowledge of the D l ^(g), D{g) can be constructed 
by the process of the addition of representations.— If there are ni basis functions of angular momentum I in the basis 
set describing Yj- a ,o/3i we have 



D(g) =®iniD l (g), 



(B8) 



where © denotes a direct sum. 

Let us now give the explicit formulas for the D l (g). In this work only s, p, and d basis functions are used, and 
hence we restrict ourselves to / = 0, 1, and 2. Since s functions just depend on the radius, ip°(f) — ip°(r), we have 



D°(g) = 1. 



(B9) 



For I = 1 there are three p functions, p\ = p x — fx(r)x, P2 = p y = fi(r)y, and p% = p z = f\(r)z, with a certain radial 
dependence fi(r)Z^ Exploiting g^ 1 = g T , we obtain p' i = Y^j=iPjSji- Thus the 3x3 representation of the rotation 
g for the p functions is 



D\ B ) 



Qxx Qxy Qxz 
Qyx Qyy Qyz 
Qzx Qzy Qzz 



(BIO) 



For I = 2 there are five d functions, where d\ = d^ z 2_ r 2 = f%{r) (3z 2 — r 2 ) / (2-\/3), ^2 = d xz = fy(r)xz, d^ = 
d yz = jiif)yz, di = d xy = f2(r)xy, and ds — d x 2_ y 2 = f%(r) (a; 2 — y 2 ) /2 with some radial dependence /2(F)- The 
transformed d functions are given as d\ — dj \fi 2 (&)] with the 5x5 representation 



D\g) = 



V%Qzy 



Qzz. 



( {3qIz ~ l) /2 VSQxxQzz 

V^QxzQzz Qxx Qzz + Q 

V^QyzQzz Qyx Qzz + Q 

V3Qx zQyz QxxQyz ~\~ QyxQxz QxyQyz ~\~ QyyQ. 

\ \/3 {Qxz ~ Qyz) /2 Qxx Qxz Qyx Qyz Qxy Qxz 



V^QzxQzy 

zxQxz QxyQzz H~~ QxzQzy QxxQzy ~\~ QxyQzx 
Qyz QyyQzz + QyzQzy QyxQzy + QyyQzx 
Qxx Qyy ~\~ Qxy Qyx 
QyyQyz QxxQxy ~ Qyx Qyy 
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Figure 13: Necessity to impose the translational symmetry on the hopping elements Hj* e ^ e . A finite spherical cluster is 
displayed. In blue-shaded areas surface effects are important. While the overlap elements S^op" are translationally symmetric, 
this is not the case for H^^ e . In particular S^q™ = S^ 1 ^^, while H^q™ 7^ -ifl as illustrated in the plot, where 
Y - S. ll. 



4. Imposing the fee space-group 

In this section we consider how to impose the fee space-group on the parameters H^^ e , extracted from the finite 
spherical fee clusters [Fig. [2J. Assuming basis functions to be real- valued, the matrix elements of a translationally 
invariant Hamiltonian H trans are symmetric and obey the relations 

rrtrans rrtrans rrtrans rrtrans /'T310 > \ 

U ia,j/3 — n (i-j)otfiP — n Qa,-(i-j)f3 ~ n -(i-j)f3,0of \ alz > 

Owing to surface effects, the translational symmetry is not fulfilled by the parameters H^^ e , as illustrated in Fig. [T5] 
Hence, although the deviations decrease with growing radius of the spheres, the translational symmetry needs to be 
enforced in order to describe a crystal. To avoid numerical errors, we impose at the same time the point-group 
symmetry Oh although that symmetry is already present due to the shape of our clusters. Concerning the notation, 
we will call the parameters conforming to the Oh point-group, the translational symmetry, and the fee space-group 
H-fao/3' ^jao/f > an d Hjafip = HjaOBi respectively. We do not need to consider the overlap, since it depends only on 
the relative position of two atoms. 

Oh point-group symmetry With Eq. (|B7|) a Hamiltonian H®£ 0j3 conforming to the point-group symmetry can be 

constructed by averaging, for a given element of H°£ p, over all H^ 1 ^ related to it by symmetry 

H? h - = — V Y \D T (g)] H sp l ler l D(g) v0 . (B13) 

Uh geO h -i,v 

Here, g runs over all No h = 48 symmetry elements of the point group O^i 81 ' 82 

Translational symmetry Using Eq. (|B12j) . the translational symmetry can be imposed by setting 



rrtrans 
n ja,0(3 



\(*Zm+b3pZ)- ( B14 ) 



Fee space-group The combined action of the Oh point-group and the translational symmetry leads to the fee space- 

fec 
ja,O0 



group symmetry. With Eqs. (|B13|) and (|B14|) we obtain the fee space-group symmetric parameters Hj Q ,o/3 = 

H fcc 

according to the prescription 

— Y Y \\D T (g)} H sp l ler Z D(g) H + \D T (g)} a H sphe S e . D (g) ). (B15) 



H fcc _ 



Rja,O0 2N 0h 

Uh e eo h »,u 
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